
setwd("/Users/aspearot/Dropbox/MaLT 2017/Draft/submission_draft/ReStat/resubmission 2/replication package/malt/")


######## 
########  Load data again
######## 

library(AER)
library(nleqslv)
library(readstata13)
library(matrixStats)
library(estimatr)

v2<-read.csv("Data/Temp/Farmer_Replication.csv",header=TRUE)	

y<-read.csv("Data/Temp/Agrovet_Replication.csv",header=TRUE)
z<-read.csv("Data/Temp/Costs_Shares_Replication.csv",header=TRUE)	
w<-read.csv("Data/Temp/Costs_Shares_Replication.csv",header=TRUE)	

a2<-as.matrix(read.csv("Data/Temp/Adoption_Replication.csv",header=TRUE))	
f2<-as.matrix(read.csv("Data/Temp/Expenditure_Replication.csv",header=TRUE))		
id<-as.matrix(read.csv("Data/Temp/ID_Replication.csv",header=TRUE))		

av_village<-as.character(tapply(as.character(v2$village_name),paste(v2$village,v2$ward,v2$district),FUN=unique))
av_ward<-as.character(tapply(as.character(v2$ward),paste(v2$village,v2$ward,v2$district),FUN=unique))
av_district<-as.character(tapply(as.character(v2$district),paste(v2$village,v2$ward,v2$district),FUN=unique))
av_market<-as.character(tapply(as.character(v2$market),paste(v2$village,v2$ward,v2$district),FUN=unique))
av_pop_market<-as.numeric(tapply(v2$census_bothsex,paste(v2$village,v2$ward,v2$district),FUN=sum,na.rm=TRUE))
popframe<-data.frame(av_village,av_ward,av_district,av_market,av_pop_market)

#Merge populations to agrovet data

y<-merge(y,popframe,all.x=T)

v2$census_bothsex<-v2$census_bothsex/sum(v2$census_bothsex,na.rm=TRUE)

z$same<-ifelse(z$survey_region_O==z$survey_region_D&z$survey_district_O==z$survey_district_D&z$survey_ward_O==z$survey_ward_D&z$survey_village_O==z$survey_village_D,1,0)


# Calculate rural and main-road distances

z$DIST_km_direct_rural<-z$DIST_km_direct*z$rural_km_share	
z$DIST_km_direct_main<-z$DIST_km_direct*(1-z$rural_km_share)

#Ensure factors are removed from variables

v2$region<-as.character(v2$region)
v2$district<-as.character(v2$district)
v2$market<-as.character(v2$market)
v2$ward<-as.character(v2$ward)
v2$village_name<-as.character(v2$village_name)

y$av_region<-as.character(y$av_region)
y$av_district<-as.character(y$av_district)
y$av_ward<-as.character(y$av_ward)
y$av_village<-as.character(y$av_village)

for(i in 1:10){
  z[,i]<-as.character(z[,i])
}

###  Construct bilateral matrix of distance.  J (rows) will be agrovets, and I (columns) will be villages

J<-nrow(y)
I<-nrow(v2)

d<-matrix(nrow=J,ncol=I) # distance matrix
r<-matrix(nrow=J,ncol=I) # rural share matrix


#construct distance matrices

for(i in 1:I){
  
  sz<-subset(z,survey_region_O==v2$region[i]&survey_district_O==v2$district[i]&survey_ward_O==v2$ward[i]&survey_village_O==v2$village[i])
  
  for(j in 1:J){	
    
    ssz<-subset(sz,survey_region_D==toupper(y$av_region[j])&survey_district_D==y$av_district[j]&survey_ward_D==y$av_ward[j]&survey_village_D==y$av_village[j])
    
    if(nrow(ssz)>0){
      d[j,i]<-ssz$DIST_km_direct[1]
      r[j,i]<-ssz$rural_km_share[1]				
    }
    if(nrow(ssz)==0){
      d[j,i]<-NA
      r[j,i]<-NA				
    }	
  }
  if(i%%100==0){print(i)}
}

minmat<-matrix(ncol=1,nrow=ncol(d))

for(i in 1:ncol(d)){
  minmat[i,1]<-min(d[,i],na.rm=TRUE)
}

rowrem<-grep(1,ifelse(rowSums(is.na(d))>500,1,0))  #agrovets with incomplete village linkages - 500 is just an arbitrarily large number chosen based on inspection of the data
colrem<-grep(1,ifelse(colSums(is.na(d))>300,1,0))  #villages with incomplete agrovet linkages - 300 is just an arbitrarily large number chosen based on inspection of the data

d2<-d[-rowrem,-colrem]  	#remove incomplete villages and agrovets
r2<-r[-rowrem,-colrem]

# farmer level adoption, expenditures, and IDs for counterfactuals

a2<-a2[,-colrem]
f2<-f2[,-colrem]
id<-id[,-colrem]

adopt_mat<-a2
omega_mat<-f2  #assigned initially as expenditures, but make into expenditure shares

for(i in 1:nrow(omega_mat)){
  omega_mat[i,]<-omega_mat[i,]/colSums(f2,na.rm=TRUE)
}

#Caclulate aggregate expenditures at the village level by taking the average within the sample and multiplying by village size

E_fert<-v2$expense_avg*v2$census_bothsex

#Agrovet fertilizer revenues

V_fert<-y$av_fert_rev_USD_w1

#define revenues as row vectors, and remove the agrovets for which there is no distance data for their village

V_fert<-as.matrix(V_fert[-rowrem])

#define adoption and expenditures as column vector.  Adoption of fertilizer includes with and without seeds.

E_fert<-t(as.matrix(E_fert[-colrem]))
E_fert<-ifelse(is.na(E_fert),0,E_fert)

####################            ######################## 
#### Calculate estimated trade costs for fertilizer ####
####################            ########################  

#Load data from MNL
xtau<-read.dta13("Data/buying_fert_MNL_binned_replication_edit.dta")
#xtau<-read.dta13("/Users/aspearot/Dropbox/MaLT 2017/Data/MaLT Project Data/Final Analysis Datasets/buying_fert_MNL_binned_032822.dta")

#calculate matrices of main road shares and rural road shares
m2<-d2*(1-r2)
s2<-d2*(r2)

#calculate elasticity adjusted trade costs estimates
mtau<-ifelse(m2,NA,NA)
mtau<-ifelse(m2>75,exp(as.numeric(xtau[9])),mtau)
mtau<-ifelse(m2>50&m2<=75,exp(as.numeric(xtau[8])),mtau)
mtau<-ifelse(m2>40&m2<=50,exp(as.numeric(xtau[7])),mtau)
mtau<-ifelse(m2>30&m2<=40,exp(as.numeric(xtau[6])),mtau)
mtau<-ifelse(m2>20&m2<=30,exp(as.numeric(xtau[5])),mtau)
mtau<-ifelse(m2>15&m2<=20,exp(as.numeric(xtau[4])),mtau)
mtau<-ifelse(m2>10&m2<=15,exp(as.numeric(xtau[3])),mtau)
mtau<-ifelse(m2>5&m2<=10,exp(as.numeric(xtau[2])),mtau)
mtau<-ifelse(m2>0&m2<=5,exp(as.numeric(xtau[1])),mtau)
mtau<-ifelse(m2<=0,1,mtau)

rtau<-ifelse(s2,NA,NA)
rtau<-ifelse(s2>20,exp(as.numeric(xtau[14])),rtau)
rtau<-ifelse(s2>15&s2<=20,exp(as.numeric(xtau[13])),rtau)
rtau<-ifelse(s2>10&s2<=15,exp(as.numeric(xtau[12])),rtau)
rtau<-ifelse(s2>5&s2<=10,exp(as.numeric(xtau[11])),rtau)
rtau<-ifelse(s2>0&s2<=5,exp(as.numeric(xtau[10])),rtau)
rtau<-ifelse(s2<=0,1,rtau)	

tau<-mtau*rtau

#### Calibration Step 1 - Normalize Total fertilizer expenditures and total fertilizer revenues to one.

E_fert2<-E_fert/sum(E_fert,na.rm=TRUE)
V_fert2<-V_fert/sum(V_fert,na.rm=TRUE)

#Define matrix for use in the contraction mapping

Ebig<-matrix(NA,nrow=nrow(d2),ncol=ncol(d2))

for(j in 1:nrow(Ebig)){
  Ebig[j,]<-E_fert2
}

Dpos<-as.matrix(ifelse(c(V_fert2)==0,0,1))

contract_start<-function(d){
  d2<-Dpos*as.matrix(d)
  dmatpre1<-d2
  
  dmat1<-matrix(dmatpre1,nrow=nrow(tau),ncol=ncol(tau))	
  #define lambda for agrovet j but without agrovet j in the numerator
  mshares1<-t(t((tau))/colSums(dmat1*(tau)))
  
  #predict the delta value for agrovet j
  dnew1<-(V_fert2/(rowSums(mshares1*Ebig)))
  
  #apply the normalization of deltas	
  dnew1<-dnew1/sum(dnew1)
  
  return(dnew1)
}

init_d0<-V_fert2
init_d<-V_fert2

for(i in 1:5000){
  
  new_d<-contract_start(init_d)
  if(i%%1==0){print(c(i,round(max(abs(new_d-init_d),na.rm=TRUE),digits=3)*1000000,cor(init_d,init_d0),max(init_d)))}
  if(max(abs(new_d-init_d),na.rm=TRUE)<1e-16){
    print(c(i,round(max(abs(new_d-init_d),na.rm=TRUE),digits=3)*1000000,cor(init_d,init_d0),max(init_d)))
    break}
  init_d<-new_d
}

share_func<-function(d){
  #d<-init_d1
  d2<-as.matrix(d)
  dmatpre1<-d2
  
  dmatpre1<-dmatpre1/sum(dmatpre1)
  
  dmat1<-matrix(dmatpre1,nrow=nrow(tau),ncol=ncol(tau))
  
  m1<-dmat1*(tau)
  mshares1<-t(t(m1)/colSums(m1))
  
  X1<-(V_fert2-rowSums(mshares1*Ebig))
  
  #SSR<-sum(X3^2,na.rm=TRUE)
  
  return(X1)		
}

#confirm solution by showing that the distance between revenues and expected expenditures at each agrovet is very small at the solution

max(abs(share_func(init_d)))

dmatpre1<-as.matrix(init_d)
dmatpre1<-dmatpre1/sum(dmatpre1)
dmat1<-matrix(dmatpre1,nrow=nrow(tau),ncol=ncol(tau))

m1<-dmat1*(tau)
mshares1<-t(t(m1)/colSums(m1))

Phi_f<-colSums(m1)

v3<-v2[-colrem,]
v3$farmers<-16

v3$Phi_f<-as.numeric(Phi_f)
v3$google_vil_city_km_std<-(v3$google_vil_city_km-mean(v3$google_vil_city_km,na.rm=TRUE))/sqrt(var(v3$google_vil_city_km,na.rm=TRUE))

h<-read.dta13("Data/Village Distances.dta")

h$market<-NULL

v3$ord<-seq(1,nrow(v3),1)
v4<-merge(v3,h,all.x=T)
v4<-v4[order(v4$ord),]

y2<-y[-rowrem,]

y2$deltas<-as.numeric(init_d)	#Quality adjusted prices, "delta", from the paper

#### RUN IV REGRESSIONS

elasreg0<-lm(log(deltas)~log(av_urea_rprice_USD_w1)+log(exper)+av_district,y2)		
elas_reduced<-(-as.numeric(coef(elasreg0)[2]))
elas_reduced

elasreg4<-ivreg(log(deltas)~log(av_urea_rprice_USD_w1)+log(exper)+av_district|log(av_urea_wprice_USD_w1)+log(exper)+av_district,data=y2)
elas_old<-(-as.numeric(coef(elasreg4)[2]))
elas_old

elasreg_IV<-ivreg(log(deltas)~log(av_urea_rprice_USD_w1)+log(exper)+av_district|log(av_urea_wprice_15_USD_w1)+log(exper)+av_district,data=subset(y2,av_pop_market>0))
elas_IV<-(-as.numeric(coef(elasreg_IV)[2]))
elas_IV

elasreg_IV_old<-ivreg(log(deltas)~log(av_urea_rprice_USD_w1)+log(exper)+av_district|log(av_pop_market)+log(exper)+av_district,data=subset(y2,av_pop_market>0))
elas_IV_old<-(-as.numeric(coef(elasreg_IV_old)[2]))
elas_IV_old		

#save for future use

write.csv(y2,"Results/IVDataset_Replication_Update.csv",row.names=FALSE,na=".")

zero<-function(l){
  ifelse(is.na(l),0,l)
}

infna<-function(l){
  ifelse(l==Inf,NA,l)
}

deltas_before<-y2$deltas
price_before<-y2$av_urea_rprice_USD_w1

#define elasticity used for analysis
elas<-elas_IV

#calculate main counterfactual change in trade costs
iceberg0<-tau^(-1/elas)-1
iceberg_half<-iceberg0/2
tau_half<-(1+iceberg_half)^(-elas)

#calculate counterfactual change in trade costs by main and rural roads separately
miceberg0<-mtau^(-1/elas)-1
miceberg_half<-miceberg0/2
riceberg0<-rtau^(-1/elas)-1
riceberg_half<-riceberg0/2	
tau_mainhalf<-(1+miceberg_half)^(-elas)*(1+riceberg0)^(-elas)	
tau_ruralhalf<-(1+miceberg0)^(-elas)*(1+riceberg_half)^(-elas)		

minmat<-matrix(ncol=6,nrow=ncol(d2))

for(i in 1:ncol(d2)){
  minmat[i,1]<-which.min(d2[,i])
  minmat[i,2]<-m2[minmat[i,1],i]
  minmat[i,3]<-s2[minmat[i,1],i]
  minmat[i,4]<-iceberg0[minmat[i,1],i]    
  minmat[i,5]<-which.max(mshares1[,i])   
  minmat[i,6]<-iceberg0[minmat[i,5],i]
}

# Calculate iceberg costs to closest agrovet villages

summary(minmat[,2])
summary(minmat[,3])
summary(minmat[,4])
summary(minmat[,6])
sum(minmat[,4]==minmat[,6])/nrow(minmat)

# Calculate average and median iceberg costs

mean(colMins(iceberg0))
median(colMins(iceberg0))
mean(colMins(miceberg0))
median(colMins(miceberg0))
mean(colMins(riceberg0))
median(colMins(riceberg0))

##### Create Dataset to compare with actual choices

#for(i in 1:nrow(v4)){

#  ag_village<-y2$av_village
#  ag_ward<-y2$av_ward
#  ag_district<-y2$av_district
#  ag_region<-y2$av_region
#  subag<-data.frame(ag_village,ag_ward,ag_district,ag_region)
#  subag$village_name<-v4$village_name[i]
#  subag$ward<-v4$ward[i]	  
#  subag$district<-v4$district[i]
#  subag$region<-v4$region[i]
#  subag$lam<-mshares1[,i]
#  subag$weighted_dist_std<-v4$weighted_dist_std[i]

#  if(i==1){bigag<-subag} 
#  if(i>1){bigag<-rbind(bigag,subag)} 

#}

#write.csv(bigag,"/Users/aspearot/Dropbox/MaLT 2017/Model/To Linux/Predicted_Choices.csv",col.names=TRUE,row.names=FALSE)	

m1_half<-dmat1*(tau_half)
mshares1_half<-t(t(m1_half)/colSums(m1_half))
v4$Phi_f_half<-as.numeric(colSums(m1_half))	

m1_ruralhalf<-dmat1*(tau_ruralhalf)
mshares1_ruralhalf<-t(t(m1_ruralhalf)/colSums(m1_ruralhalf))
v4$Phi_f_ruralhalf<-as.numeric(colSums(m1_ruralhalf))

m1_mainhalf<-dmat1*(tau_mainhalf)
mshares1_mainhalf<-t(t(m1_mainhalf)/colSums(m1_mainhalf))
v4$Phi_f_mainhalf<-as.numeric(colSums(m1_mainhalf))	

sv4<-v4[,c("region","district","ward","village_name","Phi_f_mainhalf","Phi_f_ruralhalf","Phi_f_half","Phi_f","weighted_dist_std","sampled")]

sx<-read.dta13("Data/MaLT Calibration Replication.dta") 
names(sx)[1]<-"region"
sx$id<-seq(1,nrow(sx),1)

sx<-merge(sx,sv4,all.x=T)

#################################################################################
############ ITERATE THROUGH PARAMETER ESTIMATES AND COUNTERFACTUALS ############
#################################################################################	

### Calculate Mark-ups


elas0<-c(elas_IV) #three elasticity estimates
mult0<-c(2) #three ratios between elasticity estimates and frechet parameter epsillon

count<-0	

### define results matrix-size will be the number of different counterfacultuals multiplied by the number of statistics (in our case, coefficient estimates, standard errors, and pvalues)
size<-18
comps<-matrix(nrow=size*length(elas0)*length(mult0),ncol=11)	  

for(elas in elas0){
  for(mult in mult0){
    
    eps<-elas*mult  #calculate epsilon
    
    dmatpre1<-as.matrix(init_d)
    dmatpre1<-dmatpre1/sum(dmatpre1)  #normalize sum of quality adjusted prices to 1 (this is not required)
    
    #Add deltas to columns
    dmat1<-matrix(dmatpre1,nrow=nrow(tau),ncol=ncol(tau))
    
    m1<-dmat1*(tau)
    mshares1<-t(t(m1)/colSums(m1))  ##calculate lambdas
    
    ### Calculate composite shares for markup equation
    
    smatpre<-mshares1*Ebig
    smat<-smatpre/rowSums(smatpre)
    
    zi<-colSums(omega_mat*adopt_mat,na.rm=TRUE)
    
    Zmat<-matrix(NA,nrow=nrow(mshares1),ncol=ncol(mshares1))
    for(j in 1:nrow(Zmat)){Zmat[j,]<-zi}
    
    #calculate predicted mark-up
    y2$mshare<-elas-elas*(eps-1)/eps*rowSums(smat*mshares1*(Zmat))
    y2$est_markup<-1/(y2$mshare)
    
    #calculate baseline marginal cost using mark-up and prices		
    y2$mc<-y2$av_urea_rprice_USD_w1/(y2$est_markup+1)
    
    #calculate quality terms
    y2$Ts<-y2$deltas/y2$av_urea_rprice_USD_w1^(-elas)
    
    
    ones<-matrix(1,ncol=ncol(tau),nrow=1)		
    
    ## calculate counterfactual changes in trade costs
    
    iceberg0<-tau^(-1/elas)-1
    iceberg_half<-iceberg0/2
    tau_half<-(1+iceberg_half)^(-elas)
    
    miceberg0<-mtau^(-1/elas)-1
    miceberg_half<-miceberg0/2
    riceberg0<-rtau^(-1/elas)-1
    riceberg_half<-riceberg0/2	
    
    tau_mainhalf<-(1+miceberg_half)^(-elas)*(1+riceberg0)^(-elas)	
    tau_ruralhalf<-(1+miceberg0)^(-elas)*(1+riceberg_half)^(-elas)		
    
    #define function that solves agrovet pricing equation
    
    pricingfunc<-function(r){
      
      r<-ifelse(r==0,NA,r)  ## error wrapper that is not used
      
      #Calculate new delta's (quality adjusted prices) at new prices
      
      dmatpre<-y2$Ts*r^(-elas)
      dmatpre<-ifelse(is.na(dmatpre),0,dmatpre)
      
      ### Calculate new lambdas and the shock
      
      dmat<-matrix(dmatpre,nrow=nrow(newtau),ncol=ncol(newtau))
      m<-dmat*(newtau)
      mshares<-t(t(m)/colSums(m))
      
      Shock<-colSums(m)/Phi_f
      
      ### Calculate new adoption
      
      adopt_mat_new<-matrix(NA,nrow=nrow(adopt_mat),ncol=ncol(adopt_mat))
      for(k in 1:nrow(adopt_mat)){adopt_mat_new[k,]<-Shock/(Shock+(1-adopt_mat[k,])/adopt_mat[k,])}
      
      ### Calculate new expenditures and expenditure shares
      
      f2_new<-matrix(NA,nrow=nrow(f2),ncol=ncol(f2))
      for(k in 1:nrow(f2_new)){f2_new[k,]<-f2[k,]*Shock^(1/eps)*(adopt_mat_new[k,]/adopt_mat[k,])^((eps-1)/eps)}
      
      omega_mat_new<-matrix(NA,nrow=nrow(omega_mat),ncol=ncol(omega_mat))
      for(i in 1:nrow(omega_mat_new)){
        omega_mat_new[i,]<-f2_new[i,]/colSums(f2_new,na.rm=TRUE)
      } 
      
      ### Calculate composite shares for markup equation
      
      zi_new<-colSums(omega_mat_new*adopt_mat_new,na.rm=TRUE)
      
      Zmat_new<-matrix(NA,nrow=nrow(mshares1),ncol=ncol(mshares1))
      for(j in 1:nrow(Zmat_new)){Zmat_new[j,]<-zi_new}			
      
      ### Calculate new aggregate expenditures
      
      Ebig_new<-matrix(NA,nrow=nrow(Ebig),ncol=ncol(Ebig))
      
      for(j in 1:nrow(Ebig_new)){
        Ebig_new[j,]<-colMeans(f2_new,na.rm=TRUE)*v3$census_bothsex
      }
      
      ### Calculate mark-ups
      
      smatpre_new<-mshares*Ebig_new
      smat_new<-smatpre_new/rowSums(smatpre_new)			
      
      mshare_new<-elas-elas*(eps-1)/eps*rowSums(smat_new*mshares*(Zmat_new))
      newmarkup<-1/(mshare_new)				
      
      X1<-(r-newmc*(1+newmarkup))
      X1<-ifelse(is.na(X1),0,X1)
      return(X1)		
      
    }
    
    #this function calculates outcomes
    
    outcomes<-function(r,cost,tn){
      
      r<-ifelse(r==0,NA,r)
      
      #Calculate new delta's (quality adjusted prices) at new prices
      
      dmatpre<-y2$Ts*r^(-elas)
      dmatpre<-ifelse(is.na(dmatpre),0,dmatpre)
      
      ### Calculate new lambdas and the shock
      
      dmat<-matrix(dmatpre,nrow=nrow(newtau),ncol=ncol(newtau))
      m<-dmat*(newtau)
      mshares<-t(t(m)/colSums(m))
      
      Shock<-colSums(m)/Phi_f
      
      ### Calculate new adoption
      
      adopt_mat_new<-matrix(NA,nrow=nrow(adopt_mat),ncol=ncol(adopt_mat))
      for(k in 1:nrow(adopt_mat)){adopt_mat_new[k,]<-Shock/(Shock+(1-adopt_mat[k,])/adopt_mat[k,])}
      
      ### Calculate new expenditures
      
      f2_new<-matrix(NA,nrow=nrow(f2),ncol=ncol(f2))
      for(k in 1:nrow(f2_new)){f2_new[k,]<-f2[k,]*Shock^(1/eps)*(adopt_mat_new[k,]/adopt_mat[k,])^((eps-1)/eps)}
      
      omega_mat_new<-matrix(NA,nrow=nrow(omega_mat),ncol=ncol(omega_mat))
      for(i in 1:nrow(omega_mat_new)){
        omega_mat_new[i,]<-f2_new[i,]/colSums(f2_new,na.rm=TRUE)
      } 
      
      ### Calculate new revenues
      
      Fi<-as.matrix(colSums(f2_new,na.rm=TRUE))
      
      ### Calculate composite shares for markup equation
      
      zi_new<-colSums(omega_mat_new*adopt_mat_new,na.rm=TRUE)
      
      Zmat_new<-matrix(NA,nrow=nrow(mshares1),ncol=ncol(mshares1))
      for(j in 1:nrow(Zmat_new)){Zmat_new[j,]<-zi_new}			
      
      ### Calculate new aggregate expenditures
      
      Ebig_new<-matrix(NA,nrow=nrow(Ebig),ncol=ncol(Ebig))
      
      for(j in 1:nrow(Ebig_new)){
        Ebig_new[j,]<-colMeans(f2_new,na.rm=TRUE)*v3$census_bothsex
      }
      
      ### Calculate mark-ups, profits
      
      smatpre_new<-mshares*Ebig_new
      Rj<-rowSums(smatpre_new)			
      smat_new<-smatpre_new/rowSums(smatpre_new)			
      
      mshare_new<-elas-elas*(eps-1)/eps*rowSums(smat_new*mshares*(Zmat_new))
      marks<-1/(mshare_new)			
      
      newprofits<-Rj*(marks/(marks+1))
      
      return(list(adopt=(adopt_mat_new),expenditures=(f2_new),markup=as.numeric(marks),sales=as.numeric(Rj),profits=as.numeric(newprofits),Phihat=as.numeric(Shock)))
      
    }
    
    # baseline equilibrium (if this does not converge immediately, there is an error)
    
    init_r<-ifelse(is.na(y2$av_urea_rprice_USD_w1),0,y2$av_urea_rprice_USD_w1)
    newmc<-y2$mc
    newtau<-tau
    ans<-nleqslv(init_r,fn=pricingfunc,control=list(allowSingular=TRUE,trace=1,maxit=10000))	
    adopt_mat_base<-outcomes(ans$x,y2$mc,newtau)$adopt
    exp_mat_base<-outcomes(ans$x,y2$mc,newtau)$expenditures
    rev_base<-outcomes(ans$x,y2$mc,newtau)$sales	
    profit_base<-outcomes(ans$x,y2$mc,newtau)$profits		
    markup_base<-outcomes(ans$x,y2$mc,newtau)$markup
    prices_base<-ans$x
    shock_base<-outcomes(ans$x,y2$mc,newtau)$Phihat
    
    # baseline counterfactual  	  
    
    init_r<-ifelse(is.na(y2$av_urea_rprice_USD_w1),0,y2$av_urea_rprice_USD_w1)
    newmc<-y2$mc
    newtau<-tau_half
    ans<-nleqslv(init_r,fn=pricingfunc,control=list(allowSingular=TRUE,trace=1,maxit=10000))	
    adopt_mat_half2<-outcomes(ans$x,y2$mc,newtau)$adopt
    exp_mat_half2<-outcomes(ans$x,y2$mc,newtau)$expenditures
    rev_half<-outcomes(ans$x,y2$mc,newtau)$sales
    profit_half<-outcomes(ans$x,y2$mc,newtau)$profits
    profit_half_noentry<-profit_half
    markup_half<-outcomes(ans$x,y2$mc,newtau)$markup
    prices_half<-ans$x
    shock_half<-outcomes(ans$x,y2$mc,newtau)$Phihat
    
    # main roads counterfactual  	  
    
    newtau<-tau_mainhalf
    ans<-nleqslv(init_r,fn=pricingfunc,control=list(allowSingular=TRUE,trace=1,maxit=10000))	
    adopt_mat_mainhalf2<-outcomes(ans$x,y2$mc,newtau)$adopt
    exp_mat_mainhalf2<-outcomes(ans$x,y2$mc,newtau)$expenditures	
    rev_mainhalf<-outcomes(ans$x,y2$mc,newtau)$sales
    markup_mainhalf<-outcomes(ans$x,y2$mc,newtau)$markup
    prices_mainhalf<-ans$x
    shock_mainhalf<-outcomes(ans$x,y2$mc,newtau)$Phihat
    shock_mainhalf_trade<-colSums(mshares1*newtau/tau)
    
    # rural roads counterfactual  
    
    newtau<-tau_ruralhalf
    ans<-nleqslv(init_r,fn=pricingfunc,control=list(allowSingular=TRUE,trace=1,maxit=10000))	
    adopt_mat_ruralhalf2<-outcomes(ans$x,y2$mc,newtau)$adopt
    exp_mat_ruralhalf2<-outcomes(ans$x,y2$mc,newtau)$expenditures				
    rev_ruralhalf<-outcomes(ans$x,y2$mc,newtau)$sales
    markup_ruralhalf<-outcomes(ans$x,y2$mc,newtau)$markup
    prices_ruralhalf<-ans$x
    shock_ruralhalf<-outcomes(ans$x,y2$mc,newtau)$Phihat		
    
    # distribution counterfactual  
    
    p1<-quantile(y2$av_trans_cost_per50kg,prob=0.01,na.rm=TRUE)
    p99<-quantile(y2$av_trans_cost_per50kg,prob=0.99,na.rm=TRUE)
    y2$av_trans_cost_per50kg_w1<-y2$av_trans_cost_per50kg
    summary(y2$av_trans_cost_per50kg_w1)
    y2$av_trans_cost_per50kg_w1<-ifelse(y2$av_trans_cost_per50kg_w1>p99,p99,y2$av_trans_cost_per50kg_w1)
    y2$av_trans_cost_per50kg_w1<-ifelse(y2$av_trans_cost_per50kg_w1<p1,p1,y2$av_trans_cost_per50kg_w1)
    summary(y2$av_trans_cost_per50kg_w1)
    
    trans_glm<-glm(av_trans_cost_per50kg_w1~log(google_vil_city_km),y2,family=poisson(link="log"))
    summary(trans_glm)
    y2$impute_tradecost<-ifelse(is.na(y2$av_trans_cost_per50kg_w1),as.numeric(predict(trans_glm,newdata=y2,type="response")),y2$av_trans_cost_per50kg_w1)/2300
    
    newmc<-y2$mc
    newmc<-y2$mc-y2$impute_tradecost/2
    newtau<-tau
    ans<-nleqslv(init_r,fn=pricingfunc,control=list(allowSingular=TRUE,trace=1,maxit=10000))
    adopt_halftrans<-outcomes(ans$x,newmc,newtau)$adopt		
    exp_halftrans<-outcomes(ans$x,newmc,newtau)$expenditures
    markup_halftrans<-outcomes(ans$x,y2$mc,newtau)$markup
    prices_halftrans<-ans$x
    shock_halftrans<-outcomes(ans$x,y2$mc,newtau)$Phihat		
    
    # Subsidy counterfactual
    
    newmc<-y2$mc*0.90
    newtau<-tau
    ans<-nleqslv(init_r,fn=pricingfunc,control=list(allowSingular=TRUE,trace=1,maxit=10000))
    adopt_fisp<-outcomes(ans$x,newmc,newtau)$adopt		
    exp_fisp<-outcomes(ans$x,newmc,newtau)$expenditures
    markup_fisp<-outcomes(ans$x,newmc,newtau)$markup		
    prices_fisp<-ans$x
    shock_fisp<-outcomes(ans$x,newmc,newtau)$Phihat			
    
    
    # Create dataset of farmer-level predictions from baseline, and allcounterfactual		
    
    N<-50
    
    res<-matrix(nrow=N*nrow(v4),ncol=21)
    
    for(i in 1:nrow(v4)){
      
      res[(1+(i-1)*N):(i*N),1]<-id[,i]
      res[(1+(i-1)*N):(i*N),2]<-adopt_mat[,i]
      res[(1+(i-1)*N):(i*N),3]<-f2[,i]
      #res[(1+(i-1)*N):(i*N),3]<-v4$weighted_dist_std[i]
      #res[(1+(i-1)*N):(i*N),4]<-v4$sampled[i]	  
      res[(1+(i-1)*N):(i*N),4]<-adopt_mat_half2[,i]
      res[(1+(i-1)*N):(i*N),5]<-exp_mat_half2[,i]	  
      res[(1+(i-1)*N):(i*N),6]<-adopt_mat_mainhalf2[,i]
      res[(1+(i-1)*N):(i*N),7]<-exp_mat_mainhalf2[,i]	  
      res[(1+(i-1)*N):(i*N),8]<-adopt_mat_ruralhalf2[,i]
      res[(1+(i-1)*N):(i*N),9]<-exp_mat_ruralhalf2[,i]	  
      res[(1+(i-1)*N):(i*N),10]<-adopt_halftrans[,i]
      res[(1+(i-1)*N):(i*N),11]<-exp_halftrans[,i]
      res[(1+(i-1)*N):(i*N),12]<-adopt_fisp[,i]
      res[(1+(i-1)*N):(i*N),13]<-exp_fisp[,i]	
      res[(1+(i-1)*N):(i*N),14]<-i
      res[(1+(i-1)*N):(i*N),15]<-shock_base[i]
      res[(1+(i-1)*N):(i*N),16]<-shock_half[i]
      res[(1+(i-1)*N):(i*N),17]<-shock_mainhalf[i]
      res[(1+(i-1)*N):(i*N),18]<-shock_ruralhalf[i]
      res[(1+(i-1)*N):(i*N),19]<-shock_halftrans[i]
      res[(1+(i-1)*N):(i*N),20]<-shock_fisp[i]
      res[(1+(i-1)*N):(i*N),21]<-v4$weighted_dist_std[i]
    }
    
    res<-as.data.frame(res)
    names(res)<-c("id","adopt_base","exp_base","adopt_half","exp_half","adopt_halfmain","exp_halfmain","adopt_halfrural","exp_halfrural","adopt_halftrans","exp_halftrans","adopt_fisp","exp_fisp","village","shock_base","shock_half","shock_mainhalf","shock_ruralhalf","shock_halftrans","shock_fisp","weighted_dist_std")
    
    sx2<-subset(res)
    
    # Run regressions of baseline, and change in log outcomes, on standardized distance
    
    base<-summary(lm_robust(log(adopt_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))		
    half<-summary(lm_robust(log(adopt_half/adopt_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    halfmain<-summary(lm_robust(log(adopt_halfmain/adopt_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    halfrural<-summary(lm_robust(log(adopt_halfrural/adopt_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    halftrans<-summary(lm_robust(log(adopt_halftrans/adopt_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    fisp<-summary(lm_robust(log(adopt_fisp/adopt_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    
    shockreg_base<-summary(lm_robust(log(shock_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))		
    shockreg_half<-summary(lm_robust(log(shock_half)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    shockreg_mainhalf<-summary(lm_robust(log(shock_mainhalf)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    shockreg_ruralhalf<-summary(lm_robust(log(shock_ruralhalf)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    shockreg_halftrans<-summary(lm_robust(log(shock_halftrans)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    shockreg_fisp<-summary(lm_robust(log(shock_fisp)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    
    baseexp<-summary(lm_robust(log(exp_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))		
    halfexp<-summary(lm_robust(log(exp_half/exp_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    halfmainexp<-summary(lm_robust(log(exp_halfmain/exp_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    halfruralexp<-summary(lm_robust(log(exp_halfrural/exp_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))		
    halftransexp<-summary(lm_robust(log(exp_halftrans/exp_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))		
    fispexp<-summary(lm_robust(log(exp_fisp/exp_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    
    shockreg2_base<-summary(lm_robust(I(-log(1+adopt_base*(shock_base-1)))~weighted_dist_std,cluster=village,sx2,se_type="stata"))		
    shockreg2_half<-summary(lm_robust(I(-log(1+adopt_base*(shock_half-1)))~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    shockreg2_mainhalf<-summary(lm_robust(I(-log(1+adopt_base*(shock_mainhalf-1)))~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    shockreg2_ruralhalf<-summary(lm_robust(I(-log(1+adopt_base*(shock_ruralhalf-1)))~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    shockreg2_halftrans<-summary(lm_robust(I(-log(1+adopt_base*(shock_halftrans-1)))~weighted_dist_std,cluster=village,sx2,se_type="stata"))
    shockreg2_fisp<-summary(lm_robust(I(-log(1+adopt_base*(shock_fisp-1)))~weighted_dist_std,cluster=village,sx2,se_type="stata"))		
    
    ## Assign to results file
    
    comps[(size*count+1):(size*count+size),1]<-elas
    comps[(size*count+1):(size*count+size),2]<-eps	
    
    comps[(size*count+1),3]<-round(base$coefficients[1,1],digits=4)
    comps[(size*count+2),3]<-round(base$coefficients[1,2],digits=4)
    comps[(size*count+3),3]<-round(base$coefficients[1,4],digits=4)
    comps[(size*count+4),3]<-round(half$coefficients[1,1],digits=4)
    comps[(size*count+5),3]<-round(half$coefficients[1,2],digits=4)
    comps[(size*count+6),3]<-round(half$coefficients[1,4],digits=4)
    comps[(size*count+7),3]<-round(halfrural$coefficients[1,1],digits=4)
    comps[(size*count+8),3]<-round(halfrural$coefficients[1,2],digits=4)
    comps[(size*count+9),3]<-round(halfrural$coefficients[1,4],digits=4)
    comps[(size*count+10),3]<-round(halfmain$coefficients[1,1],digits=4)
    comps[(size*count+11),3]<-round(halfmain$coefficients[1,2],digits=4)
    comps[(size*count+12),3]<-round(halfmain$coefficients[1,4],digits=4)	
    comps[(size*count+13),3]<-round(halftrans$coefficients[1,1],digits=4)
    comps[(size*count+14),3]<-round(halftrans$coefficients[1,2],digits=4)
    comps[(size*count+15),3]<-round(halftrans$coefficients[1,4],digits=4)	
    comps[(size*count+16),3]<-round(fisp$coefficients[1,1],digits=4)
    comps[(size*count+17),3]<-round(fisp$coefficients[1,2],digits=4)
    comps[(size*count+18),3]<-round(fisp$coefficients[1,4],digits=4)	
    
    comps[(size*count+1),4]<-round(base$coefficients[2,1],digits=4)
    comps[(size*count+2),4]<-round(base$coefficients[2,2],digits=4)
    comps[(size*count+3),4]<-round(base$coefficients[2,4],digits=4)
    comps[(size*count+4),4]<-round(half$coefficients[2,1],digits=4)
    comps[(size*count+5),4]<-round(half$coefficients[2,2],digits=4)
    comps[(size*count+6),4]<-round(half$coefficients[2,4],digits=4)
    comps[(size*count+7),4]<-round(halfrural$coefficients[2,1],digits=4)
    comps[(size*count+8),4]<-round(halfrural$coefficients[2,2],digits=4)
    comps[(size*count+9),4]<-round(halfrural$coefficients[2,4],digits=4)
    comps[(size*count+10),4]<-round(halfmain$coefficients[2,1],digits=4)
    comps[(size*count+11),4]<-round(halfmain$coefficients[2,2],digits=4)
    comps[(size*count+12),4]<-round(halfmain$coefficients[2,4],digits=4)	
    comps[(size*count+13),4]<-round(halftrans$coefficients[2,1],digits=4)
    comps[(size*count+14),4]<-round(halftrans$coefficients[2,2],digits=4)
    comps[(size*count+15),4]<-round(halftrans$coefficients[2,4],digits=4)	
    comps[(size*count+16),4]<-round(fisp$coefficients[2,1],digits=4)
    comps[(size*count+17),4]<-round(fisp$coefficients[2,2],digits=4)
    comps[(size*count+18),4]<-round(fisp$coefficients[2,4],digits=4)	
    
    comps[(size*count+1),5]<-round(baseexp$coefficients[1,1],digits=4)
    comps[(size*count+2),5]<-round(baseexp$coefficients[1,2],digits=4)
    comps[(size*count+3),5]<-round(baseexp$coefficients[1,4],digits=4)
    comps[(size*count+4),5]<-round(halfexp$coefficients[1,1],digits=4)
    comps[(size*count+5),5]<-round(halfexp$coefficients[1,2],digits=4)
    comps[(size*count+6),5]<-round(halfexp$coefficients[1,4],digits=4)
    comps[(size*count+7),5]<-round(halfruralexp$coefficients[1,1],digits=4)
    comps[(size*count+8),5]<-round(halfruralexp$coefficients[1,2],digits=4)
    comps[(size*count+9),5]<-round(halfruralexp$coefficients[1,4],digits=4)
    comps[(size*count+10),5]<-round(halfmainexp$coefficients[1,1],digits=4)
    comps[(size*count+11),5]<-round(halfmainexp$coefficients[1,2],digits=4)
    comps[(size*count+12),5]<-round(halfmainexp$coefficients[1,4],digits=4)	
    comps[(size*count+13),5]<-round(halftransexp$coefficients[1,1],digits=4)
    comps[(size*count+14),5]<-round(halftransexp$coefficients[1,2],digits=4)
    comps[(size*count+15),5]<-round(halftransexp$coefficients[1,4],digits=4)	
    comps[(size*count+16),5]<-round(fispexp$coefficients[1,1],digits=4)
    comps[(size*count+17),5]<-round(fispexp$coefficients[1,2],digits=4)
    comps[(size*count+18),5]<-round(fispexp$coefficients[1,4],digits=4)	
    
    comps[(size*count+1),6]<-round(baseexp$coefficients[2,1],digits=4)
    comps[(size*count+2),6]<-round(baseexp$coefficients[2,2],digits=4)
    comps[(size*count+3),6]<-round(baseexp$coefficients[2,4],digits=4)
    comps[(size*count+4),6]<-round(halfexp$coefficients[2,1],digits=4)
    comps[(size*count+5),6]<-round(halfexp$coefficients[2,2],digits=4)
    comps[(size*count+6),6]<-round(halfexp$coefficients[2,4],digits=4)
    comps[(size*count+7),6]<-round(halfruralexp$coefficients[2,1],digits=4)
    comps[(size*count+8),6]<-round(halfruralexp$coefficients[2,2],digits=4)
    comps[(size*count+9),6]<-round(halfruralexp$coefficients[2,4],digits=4)
    comps[(size*count+10),6]<-round(halfmainexp$coefficients[2,1],digits=4)
    comps[(size*count+11),6]<-round(halfmainexp$coefficients[2,2],digits=4)
    comps[(size*count+12),6]<-round(halfmainexp$coefficients[2,4],digits=4)	
    comps[(size*count+13),6]<-round(halftransexp$coefficients[2,1],digits=4)
    comps[(size*count+14),6]<-round(halftransexp$coefficients[2,2],digits=4)
    comps[(size*count+15),6]<-round(halftransexp$coefficients[2,4],digits=4)	
    comps[(size*count+16),6]<-round(fispexp$coefficients[2,1],digits=4)
    comps[(size*count+17),6]<-round(fispexp$coefficients[2,2],digits=4)
    comps[(size*count+18),6]<-round(fispexp$coefficients[2,4],digits=4)	
    
    comps[(size*count+1),7]<-round(shockreg_base$coefficients[1,1],digits=4)
    comps[(size*count+2),7]<-round(shockreg_base$coefficients[1,2],digits=4)
    comps[(size*count+3),7]<-round(shockreg_base$coefficients[1,4],digits=4)
    comps[(size*count+4),7]<-round(shockreg_half$coefficients[1,1],digits=4)
    comps[(size*count+5),7]<-round(shockreg_half$coefficients[1,2],digits=4)
    comps[(size*count+6),7]<-round(shockreg_half$coefficients[1,4],digits=4)
    comps[(size*count+7),7]<-round(shockreg_ruralhalf$coefficients[1,1],digits=4)
    comps[(size*count+8),7]<-round(shockreg_ruralhalf$coefficients[1,2],digits=4)
    comps[(size*count+9),7]<-round(shockreg_ruralhalf$coefficients[1,4],digits=4)
    comps[(size*count+10),7]<-round(shockreg_mainhalf$coefficients[1,1],digits=4)
    comps[(size*count+11),7]<-round(shockreg_mainhalf$coefficients[1,2],digits=4)
    comps[(size*count+12),7]<-round(shockreg_mainhalf$coefficients[1,4],digits=4)	
    comps[(size*count+13),7]<-round(shockreg_halftrans$coefficients[1,1],digits=4)
    comps[(size*count+14),7]<-round(shockreg_halftrans$coefficients[1,2],digits=4)
    comps[(size*count+15),7]<-round(shockreg_halftrans$coefficients[1,4],digits=4)	
    comps[(size*count+16),7]<-round(shockreg_fisp$coefficients[1,1],digits=4)
    comps[(size*count+17),7]<-round(shockreg_fisp$coefficients[1,2],digits=4)
    comps[(size*count+18),7]<-round(shockreg_fisp$coefficients[1,4],digits=4)	
    
    comps[(size*count+1),8]<-round(shockreg_base$coefficients[2,1],digits=4)
    comps[(size*count+2),8]<-round(shockreg_base$coefficients[2,2],digits=4)
    comps[(size*count+3),8]<-round(shockreg_base$coefficients[2,4],digits=4)
    comps[(size*count+4),8]<-round(shockreg_half$coefficients[2,1],digits=4)
    comps[(size*count+5),8]<-round(shockreg_half$coefficients[2,2],digits=4)
    comps[(size*count+6),8]<-round(shockreg_half$coefficients[2,4],digits=4)
    comps[(size*count+7),8]<-round(shockreg_ruralhalf$coefficients[2,1],digits=4)
    comps[(size*count+8),8]<-round(shockreg_ruralhalf$coefficients[2,2],digits=4)
    comps[(size*count+9),8]<-round(shockreg_ruralhalf$coefficients[2,4],digits=4)
    comps[(size*count+10),8]<-round(shockreg_mainhalf$coefficients[2,1],digits=4)
    comps[(size*count+11),8]<-round(shockreg_mainhalf$coefficients[2,2],digits=4)
    comps[(size*count+12),8]<-round(shockreg_mainhalf$coefficients[2,4],digits=4)	
    comps[(size*count+13),8]<-round(shockreg_halftrans$coefficients[2,1],digits=4)
    comps[(size*count+14),8]<-round(shockreg_halftrans$coefficients[2,2],digits=4)
    comps[(size*count+15),8]<-round(shockreg_halftrans$coefficients[2,4],digits=4)	
    comps[(size*count+16),8]<-round(shockreg_fisp$coefficients[2,1],digits=4)
    comps[(size*count+17),8]<-round(shockreg_fisp$coefficients[2,2],digits=4)
    comps[(size*count+18),8]<-round(shockreg_fisp$coefficients[2,4],digits=4)			
    
    comps[(size*count+1),9]<-round(shockreg2_base$coefficients[1,1],digits=4)
    comps[(size*count+2),9]<-round(shockreg2_base$coefficients[1,2],digits=4)
    comps[(size*count+3),9]<-round(shockreg2_base$coefficients[1,4],digits=4)
    comps[(size*count+4),9]<-round(shockreg2_half$coefficients[1,1],digits=4)
    comps[(size*count+5),9]<-round(shockreg2_half$coefficients[1,2],digits=4)
    comps[(size*count+6),9]<-round(shockreg2_half$coefficients[1,4],digits=4)
    comps[(size*count+7),9]<-round(shockreg2_ruralhalf$coefficients[1,1],digits=4)
    comps[(size*count+8),9]<-round(shockreg2_ruralhalf$coefficients[1,2],digits=4)
    comps[(size*count+9),9]<-round(shockreg2_ruralhalf$coefficients[1,4],digits=4)
    comps[(size*count+10),9]<-round(shockreg2_mainhalf$coefficients[1,1],digits=4)
    comps[(size*count+11),9]<-round(shockreg2_mainhalf$coefficients[1,2],digits=4)
    comps[(size*count+12),9]<-round(shockreg2_mainhalf$coefficients[1,4],digits=4)	
    comps[(size*count+13),9]<-round(shockreg2_halftrans$coefficients[1,1],digits=4)
    comps[(size*count+14),9]<-round(shockreg2_halftrans$coefficients[1,2],digits=4)
    comps[(size*count+15),9]<-round(shockreg2_halftrans$coefficients[1,4],digits=4)	
    comps[(size*count+16),9]<-round(shockreg2_fisp$coefficients[1,1],digits=4)
    comps[(size*count+17),9]<-round(shockreg2_fisp$coefficients[1,2],digits=4)
    comps[(size*count+18),9]<-round(shockreg2_fisp$coefficients[1,4],digits=4)	
    
    comps[(size*count+1),10]<-round(shockreg2_base$coefficients[2,1],digits=4)
    comps[(size*count+2),10]<-round(shockreg2_base$coefficients[2,2],digits=4)
    comps[(size*count+3),10]<-round(shockreg2_base$coefficients[2,4],digits=4)
    comps[(size*count+4),10]<-round(shockreg2_half$coefficients[2,1],digits=4)
    comps[(size*count+5),10]<-round(shockreg2_half$coefficients[2,2],digits=4)
    comps[(size*count+6),10]<-round(shockreg2_half$coefficients[2,4],digits=4)
    comps[(size*count+7),10]<-round(shockreg2_ruralhalf$coefficients[2,1],digits=4)
    comps[(size*count+8),10]<-round(shockreg2_ruralhalf$coefficients[2,2],digits=4)
    comps[(size*count+9),10]<-round(shockreg2_ruralhalf$coefficients[2,4],digits=4)
    comps[(size*count+10),10]<-round(shockreg2_mainhalf$coefficients[2,1],digits=4)
    comps[(size*count+11),10]<-round(shockreg2_mainhalf$coefficients[2,2],digits=4)
    comps[(size*count+12),10]<-round(shockreg2_mainhalf$coefficients[2,4],digits=4)	
    comps[(size*count+13),10]<-round(shockreg2_halftrans$coefficients[2,1],digits=4)
    comps[(size*count+14),10]<-round(shockreg2_halftrans$coefficients[2,2],digits=4)
    comps[(size*count+15),10]<-round(shockreg2_halftrans$coefficients[2,4],digits=4)	
    comps[(size*count+16),10]<-round(shockreg2_fisp$coefficients[2,1],digits=4)
    comps[(size*count+17),10]<-round(shockreg2_fisp$coefficients[2,2],digits=4)
    comps[(size*count+18),10]<-round(shockreg2_fisp$coefficients[2,4],digits=4)			
    
    
    comps[(size*count+1),11]<-"base"
    comps[(size*count+2),11]<-"base"
    comps[(size*count+3),11]<-"base"
    comps[(size*count+4),11]<-"half"
    comps[(size*count+5),11]<-"half"
    comps[(size*count+6),11]<-"half"
    comps[(size*count+7),11]<-"rural"
    comps[(size*count+8),11]<-"rural"	
    comps[(size*count+9),11]<-"rural"
    comps[(size*count+10),11]<-"main"
    comps[(size*count+11),11]<-"main"	
    comps[(size*count+12),11]<-"main"	
    comps[(size*count+13),11]<-"dist"
    comps[(size*count+14),11]<-"dist"	
    comps[(size*count+15),11]<-"dist"
    comps[(size*count+16),11]<-"fisp"
    comps[(size*count+17),11]<-"fisp"	
    comps[(size*count+18),11]<-"fisp"	
    
    count<-count+1
    
  }}



#calculate profits, normalized and raw
y2$profits<-V_fert2*y2$est_markup/(y2$est_markup+1)	
y2$profits_raw<-V_fert*y2$est_markup/(y2$est_markup+1)	

#save profits from baseline equilibrium, and the equilibrium with half trade costs
y2$profits_base<-profit_base
y2$profits_half<-profit_half
				

names(h)[1:4]<-c("av_district","av_ward","av_village","av_region")
y3<-merge(y2,h,all.x=T)

#run regression of marginal cost on remoteness, to use when adding a hypothetical agrovet in a particular location
mcreg<-lm(mc~weighted_dist_std,y3)


#### Create distance matrix to call during the loops to increase efficiency.

big_tau<-matrix(nrow=ncol(tau),ncol=ncol(tau))
big_r<-matrix(nrow=ncol(tau),ncol=ncol(tau))

for(i in 1:ncol(tau)){
  
  sz<-subset(z,survey_region_O==v3$region[i]&survey_district_O==v3$district[i]&survey_ward_O==v3$ward[i]&survey_village_O==v3$village[i])
  
  for(j in 1:ncol(tau)){	
    
    ssz<-subset(sz,survey_region_D==toupper(v3$region[j])&survey_district_D==v3$district[j]&survey_ward_D==v3$ward[j]&survey_village_D==v3$village[j])
    
    if(nrow(ssz)>0){
      big_tau[i,j]<-ssz$DIST_km_direct[1]
      big_r[i,j]<-ssz$rural_km_share[1]
    }
    if(nrow(ssz)==0){
      big_tau[i,j]<-NA
      big_r[i,j]<-NA
    }			
     
  }
  print(i)
}

#############  ITERATE THROUGH ALL VILLAGES TO ADD A HYPOTHETICAL ENTRANT  #############

for(i in 1:ncol(tau)){
  
  #define and populate distance and rural share of transport for a new entrant in village i
  
  d2X<-matrix(nrow=1,ncol=ncol(tau))
  d2X<-big_tau[i,]
  r2X<-matrix(nrow=1,ncol=ncol(tau))  
  r2X<-big_r[i,]  
  
  #Calculate the elasticity adjusted iceberg for villages to reach this new agrovet using main road distances and rural road distances
  
  m2X<-d2X*(1-r2X)
  s2X<-d2X*(r2X)

  mtauX<-ifelse(m2X,NA,NA)
  mtauX<-ifelse(m2X>75,exp(as.numeric(xtau[9])),mtauX)
  mtauX<-ifelse(m2X>50&m2X<=75,exp(as.numeric(xtau[8])),mtauX)
  mtauX<-ifelse(m2X>40&m2X<=50,exp(as.numeric(xtau[7])),mtauX)
  mtauX<-ifelse(m2X>30&m2X<=40,exp(as.numeric(xtau[6])),mtauX)
  mtauX<-ifelse(m2X>20&m2X<=30,exp(as.numeric(xtau[5])),mtauX)
  mtauX<-ifelse(m2X>15&m2X<=20,exp(as.numeric(xtau[4])),mtauX)
  mtauX<-ifelse(m2X>10&m2X<=15,exp(as.numeric(xtau[3])),mtauX)
  mtauX<-ifelse(m2X>5&m2X<=10,exp(as.numeric(xtau[2])),mtauX)
  mtauX<-ifelse(m2X>0&m2X<=5,exp(as.numeric(xtau[1])),mtauX)
  mtauX<-ifelse(m2X<=0,1,mtauX)
  
  rtauX<-ifelse(s2X,NA,NA)
  rtauX<-ifelse(s2X>20,exp(as.numeric(xtau[14])),rtauX)
  rtauX<-ifelse(s2X>15&s2X<=20,exp(as.numeric(xtau[13])),rtauX)
  rtauX<-ifelse(s2X>10&s2X<=15,exp(as.numeric(xtau[12])),rtauX)
  rtauX<-ifelse(s2X>5&s2X<=10,exp(as.numeric(xtau[11])),rtauX)
  rtauX<-ifelse(s2X>0&s2X<=5,exp(as.numeric(xtau[10])),rtauX)
  rtauX<-ifelse(s2X<=0,1,rtauX)	
  
  tauX<-mtauX*rtauX
  
  #continue if distances to this village are fully populated
  if(sum(is.na(tauX))==0){
  
    
  #iterate over four cases of observed agrovet quality - 25th, 50th, and 75th percentiles, and the mean (b=4)    
  for(b in 1:4){
    
    #quality vector
    Tentrants<-c(as.numeric(quantile(y2$Ts,prob=c(.25,0.5,0.75),na.rm=TRUE)),mean(y2$Ts,na.rm=TRUE))
    
    pricingfuncX<-function(r){
      
      r<-ifelse(r==0,NA,r)
      
      #Calculate new delta's (quality adjusted prices) at new prices, and with the vector qualities with the new entrant, TsX
      
      dmatpre<-TsX*r^(-elas)
      dmatpre<-ifelse(is.na(dmatpre),0,dmatpre)
      
      ### Calculate new lambdas and the shock relative to the baseline equilibrium
      
      dmat<-matrix(dmatpre,nrow=nrow(newtau),ncol=ncol(newtau))
      m<-dmat*(newtau)
      
      mshares<-t(t(m)/colSums(m))
      
      Shock<-colSums(m)/Phi_f
      
      ### Calculate new adoption
      
      adopt_mat_new<-matrix(NA,nrow=nrow(adopt_mat),ncol=ncol(adopt_mat))
      for(k in 1:nrow(adopt_mat)){adopt_mat_new[k,]<-Shock/(Shock+(1-adopt_mat[k,])/adopt_mat[k,])}
      
      ### Calculate new expenditures
      
      f2_new<-matrix(NA,nrow=nrow(f2),ncol=ncol(f2))
      for(k in 1:nrow(f2_new)){f2_new[k,]<-f2[k,]*Shock^(1/eps)*(adopt_mat_new[k,]/adopt_mat[k,])^((eps-1)/eps)}
      
      omega_mat_new<-matrix(NA,nrow=nrow(omega_mat),ncol=ncol(omega_mat))
      for(i in 1:nrow(omega_mat_new)){
        omega_mat_new[i,]<-f2_new[i,]/colSums(f2_new,na.rm=TRUE)
      } 
      
      ### Calculate composite shares for markup equation
      
      zi_new<-1+(eps-1)/eps-(eps-1)/eps*colSums(omega_mat_new*adopt_mat_new,na.rm=TRUE)
      
      Zmat_new<-matrix(NA,nrow=nrow(mshares),ncol=ncol(mshares))
      for(j in 1:nrow(Zmat_new)){Zmat_new[j,]<-zi_new}			
      
      ### Calculate new aggregate expenditures
      
      Ebig_new<-matrix(NA,nrow=nrow(newtau),ncol=ncol(newtau))
      
      for(j in 1:nrow(Ebig_new)){
        Ebig_new[j,]<-colMeans(f2_new,na.rm=TRUE)*v3$census_bothsex
      }
      
      smatpre_new<-mshares*Ebig_new
      smat_new<-smatpre_new/rowSums(smatpre_new)			
      
      mshare_new<-1-rowSums(smat_new*mshares*(Zmat_new-1))
      newmarkup<-1/(elas*mshare_new)				
      
      X1<-(r-newmcX*(1+newmarkup))
      X1<-ifelse(is.na(X1),0,X1)
      return(X1)		
    }		
    
    outcomesX<-function(r){
      
      #r<-ans$x
      #tn<-newtau
      
      r<-ifelse(r==0,NA,r)
      
      #Calculate new delta's (quality adjusted prices) at new prices
      
      dmatpre<-TsX*r^(-elas)
      dmatpre<-ifelse(is.na(dmatpre),0,dmatpre)
      
      ### Calculate new lambdas and the shock
      
      dmat<-matrix(dmatpre,nrow=nrow(newtau),ncol=ncol(newtau))
      m<-dmat*(newtau)
      
      mshares<-t(t(m)/colSums(m))
      
      Shock<-colSums(m)/Phi_f
      
      ### Calculate new adoption
      
      adopt_mat_new<-matrix(NA,nrow=nrow(adopt_mat),ncol=ncol(adopt_mat))
      for(k in 1:nrow(adopt_mat)){adopt_mat_new[k,]<-Shock/(Shock+(1-adopt_mat[k,])/adopt_mat[k,])}
      
      ### Calculate new expenditures
      
      f2_new<-matrix(NA,nrow=nrow(f2),ncol=ncol(f2))
      for(k in 1:nrow(f2_new)){f2_new[k,]<-f2[k,]*Shock^(1/eps)*(adopt_mat_new[k,]/adopt_mat[k,])^((eps-1)/eps)}
      
      omega_mat_new<-matrix(NA,nrow=nrow(omega_mat),ncol=ncol(omega_mat))
      for(i in 1:nrow(omega_mat_new)){
        omega_mat_new[i,]<-f2_new[i,]/colSums(f2_new,na.rm=TRUE)
      } 
      
      ### Calculate new revenues
      
      Fi<-as.matrix(colSums(f2_new,na.rm=TRUE))
      
      ### Calculate composite shares for markup equation
      
      zi_new<-1+(eps-1)/eps-(eps-1)/eps*colSums(omega_mat_new*adopt_mat_new,na.rm=TRUE)
      
      Zmat_new<-matrix(NA,nrow=nrow(mshares),ncol=ncol(mshares))
      for(j in 1:nrow(Zmat_new)){Zmat_new[j,]<-zi_new}			
      
      ### Calculate new aggregate expenditures
      
      Ebig_new<-matrix(NA,nrow=nrow(newtau),ncol=ncol(newtau))
      
      for(j in 1:nrow(Ebig_new)){
        Ebig_new[j,]<-colMeans(f2_new,na.rm=TRUE)*v3$census_bothsex
      }
      
      smatpre_new<-mshares*Ebig_new
      Rj<-rowSums(smatpre_new)			
      smat_new<-smatpre_new/rowSums(smatpre_new)			
      
      mshare_new<-1-rowSums(smat_new*mshares*(Zmat_new-1))
      marks<-1/(elas*mshare_new)				
      
      newprofits<-Rj*(marks/(marks+1))
      
      return(list(adopt=(adopt_mat_new),expenditures=(f2_new),markup=as.numeric(marks),sales=as.numeric(Rj),profits=as.numeric(newprofits)))
      
    }
    
    
    new_TsX<-Tentrants[b]   #asssign quality to new entrant
    TsX<-as.matrix(c(y2$Ts,new_TsX))   #add to vector of qualities

    new_mc<-as.numeric(coef(mcreg)[1])+v4$weighted_dist_std[i]*as.numeric(coef(mcreg)[2])  #assign new marginal cost using the mc on remoteness regression
    new_mc<-ifelse(is.na(new_mc),mean(subset(y2)$mc,na.rm=TRUE),new_mc)   #assign average if prediction is missing
    newmcX<-as.matrix(c(y2$mc,new_mc))  #add to marginal cost vector
    
    new_price<-median(subset(y2,av_region==v3$region[i])$av_urea_rprice_USD_w1,na.rm=TRUE)  #determine starting value for new price
    #new_price<-ifelse(is.na(new_price),median(subset(y2)$av_urea_rprice_USD_w1),new_price)		#error handling
    
    init_r<-c(ifelse(is.na(y2$av_urea_rprice_USD_w1),0,y2$av_urea_rprice_USD_w1),new_price)  #add new price to vector of starting values for prices
    init_r[length(init_r)]<-ifelse(is.na(init_r[length(init_r)]),median(subset(y2)$av_urea_rprice_USD_w1),init_r[length(init_r)])	  #error handling			
    
    newtau<-rbind(tau,tauX)  #combine existing trade costs to agrovets with the new hypothetical agrovet.
    
    #### Entry counterfactuals at baseline trade costs (using an iterative solution method)
    
    tol<-1e-8
    
    for(count in 1:1000){
      init_r1<-init_r-pricingfuncX(init_r)
      contract<-max(abs(init_r1-init_r))
      if(contract<tol){
        #print(contract)
        break}
      init_r<-init_r1
    }
    if(count==1000){
      init_r0<-c(ifelse(is.na(y2$av_urea_rprice_USD_w1),0,y2$av_urea_rprice_USD_w1),median(subset(y2,av_region==v3$region[i])$av_urea_rprice_USD_w1,na.rm=TRUE))
      init_r0[length(init_r0)]<-ifelse(is.na(init_r0[length(init_r0)]),median(subset(y2)$av_urea_rprice_USD_w1),init_r0[length(init_r0)])
      ans<-nleqslv(init_r0,fn=pricingfuncX,control=list(allowSingular=TRUE,trace=1,maxit=10000))		
      contract<-max(abs(ans$fvec))
      init_r<-ans$x
    }
    
    init_r0<-init_r
    outfile<-outcomesX(init_r)
    
    #cut trade costs in half, and find new equilibrium
    
    newiceberg0<-newtau^(-1/elas)-1
    newiceberg_half<-newiceberg0/2
    newtau_half<-(1+newiceberg_half)^(-elas)		
    newtau<-newtau_half
    
    for(count in 1:1000){
      init_r1<-init_r-pricingfuncX(init_r)
      contract<-max(abs(init_r1-init_r))
      if(contract<tol){
        #print(contract)
        break}
      init_r<-init_r1
    }
    if(count==1000){
      ans<-nleqslv(init_r0,fn=pricingfuncX,control=list(allowSingular=TRUE,trace=1,maxit=10000))		
      contract<-max(abs(ans$fvec))
      init_r<-ans$x
    }		
    
    outfile_half<-outcomesX(init_r)		
    
    region<-v3$region[i]
    district<-v3$district[i]
    ward<-v3$ward[i]
    village<-v3$village[i]
    
    profit_inc<-sum(outfile$profits[1:(length(outfile$profits)-1)],na.rm=TRUE)  #Profit of incumbents
    profit_half_inc<-sum(outfile_half$profits[1:(length(outfile$profits)-1)])  #Profit of incumbents with half trade costs
    profit<-outfile$profits[length(outfile$profits)]  #Profit of entrant
    profit_half<-outfile_half$profits[length(outfile_half$profits)]  #Profit of entrant with half trade costs
    
    if(i*b==1){entry<-data.frame(region,district,ward,village,profit,profit_half,profit_inc,profit_half_inc,b,Tentrants[b])}
    if(i*b>1){entry<-rbind(entry,data.frame(region,district,ward,village,profit,profit_half,profit_inc,profit_half_inc,b,Tentrants[b]))}
  
  }  
    
  }
  
  if(i%%100==0){print(paste("############# ITERATION",i,"#############"))
    
    write.csv(entry,"Results/Replication_Log_Update.csv",row.names=FALSE)
    
    }  
}






entry<-read.csv("Results/Replication_Log_Update.csv",header=TRUE)
names(entry)[4]<-"village_name"

h<-read.dta13("/Users/aspearot/Dropbox/MaLT 2017/Data/MaLT Project Data/Final Analysis Datasets/Village Distances.dta")

h$market<-NULL

v5<-merge(entry,h,all.x=T)

summary(lm(log(profit)~weighted_dist_std,subset(v5,b==4)))
summary(lm(log(profit_half)~weighted_dist_std,subset(v5,b==4)))


